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ABSTRACT 

Previous work reported a bar signature in color-selected IRAS variable 
stars. Here, we estimate the source density of these variables while consistently 
accounting for spatial incompleteness in data using a likelihood approach. The 
existence of the bar is confirmed with shoulder at a ~ 4 kpc, axis ratio a : b = 2.2 
-2.7 and position angle of 19° ± 1° degrees. The ratio of non-axisymmetric to 
axisymmetric components gives similar estimate for the bar size a = 3.3 ± 0.1 
kpc and position angle 00 = 24° ± 2°. We estimate a scale length 4.00 ± 0.55 
kpc for the IRAS variable population, suggesting that these stars represent the 
old disk population. 

We use this density reconstruction to estimate the optical depth to 
microlensing for the large-scale bar in the Galactic disk. We hnd an 
enhancement over an equivalent axisymmetric disk by 30% but still too small to 
account for the MACHO result. In addition, we predict a signihcant asymmetry 
at positive and negative longitudes along lines of sight through the end of the 
bar (|/| ^ 30°) with optical depths comparable to that in Baade’s window. An 
infrared microlensing survey may be a sensitive tool for detecting or constraining 
structural asymmetries. 

More generally, this is a pilot study for Bayesian star count analyses. 
Bayesian approach allows the assessment of prior probabilities to the unknown 
parameters of the model; the resulting likelihood function is straightforwardly 
modihed to incorporate all available data. However, this method requires the 
evaluation of multidimensional density functions over the data and optimization 
of the function over a parameter space. We address the resulting computational 
extremization problem with a hybrid use of a directed search algorithm which 
locates the global maximum and the conjugate gradient method which converges 



quickly near a likelihood maximum. Both methods are parallelizable and 
therefore of potential use with very large databases. 

Subject headings: Galaxy: structure — stars: variables: other — stars: 
AGB and post-AGB — gravitational lensing — methods: statistical 
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1. Introduction 

Weinberg (1992, Paper I) identified color-selected variables in the IRAS 
Point Source Catalog (PSC) with AGB stars based on color consistency and 
the circumstantial sensitivity of the IRAS survey to long-period variables (cf. 
Harmon & Gilmore 1988). These were then used as rough standard candles 
to infer a large-scale asymmetry in the stellar distribution. The identification 
of IRAS variables with AGB stars was strengthened by an in-depth study of 
a bright subset (Allen, Kleinmann & Weinberg 1993). Carbon-selected AGB 
stars (carbon stars) have also proven to be effective tracers (see e.g. Metzger & 
Schechter 1994). Advantages of AGB tracers are reviewed in Weinberg (1994). 
In general, standard candle analyses have the advantage over flux or star count 
analyses in providing direct information about the three-dimensional structure 
of the Galaxy. However, uncertainties in their selection and intrinsic properties 
may bias any inference and, especially for the IRAS-selected sample, the census 
is incomplete. 

Paper I described an approach to large-scale Galactic structure using a 
star count analysis which allows the information to be reconstructed and 
possibly corrected in the observer’s coordinate system before translating to a 
Galactocentric system. Unfortunately, this translation approach is only natural 
if the coverage is complete and suffered in application to the IRAS sample 
because of spatial gaps due to an incomplete second full-sky epoch. Here, we 
present the results of a different approach to the problem: the direct density 
estimation by maximum likelihood. A Bayesian density estimation has the 
advantage of directly incorporating selection effects and missing data. 

The number of ongoing surveys that bear on Galactic structure—SDSS, 
2MASS, DENIS—which at various stages will have surveyed parts of the sky is a 
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second motivation for this study; there is a need for a systematic method suited 
to inferential studies using possibly incomplete data from many wave bands. 
Recent analyses (e.g. Bahcall & Soneira 1980 in the optical; Wainscoat et ah 
1992 in the infrared) have modeled the Galactic components with standard 
prohles and structural parameters chosen to provide a match to star count data. 
To explore the structural parameters themselves, we propose a Bayesian density 
estimation technique to treat data from scattered helds during the survey and 
to easily incorporate data from wave bands. Conceptually, this approach is 
midway between a classical inversion and modeling. 

The hrst part of the paper describes and characterizes the method. More 
specihcally, §|^ reviews the IRAS selection procedure described in Paper I 
and motivates the approach. The new analysis based on statistical density 
estimation is presented in §|^ and precisely dehned in §^. The second part of the 
paper describes Monte-Carlo tests and the results of applying the method to the 
IRAS data (§^. We conclude in with a summary and discussion. 

2. IRAS source selection 

The analysis in Paper I was based on the variables selected in the IRAS 
Point Source Catalog (1988) by both color and Pvar- Following the source 
selection procedure described in Paper I, we selected stars from IRAS Point 
Source Catalog with F 12 > 2 Jy and variability flag P^ar > 98%. Although the 
flux limit reduces the confusion in source identihcation toward the center of the 
Galaxy, it also restricts the sensitivity to distant sources. The limiting distance 
to a star (d) is estimated using a simple exponential layer with vertical scale 
height h and mid-plane extinction coefficient Ki 2 '. 


m = M + 5 Igd - 5 + ATia h (1 - / sin |6|. 


( 1 ) 
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For a typical AGB star {L = 3500Lq, see Appendix A) and K 12 = 0.18 
kpc“^, the limiting distance in the plane is Rum = 7 kpc. We assume that the 
extinction is dominated by the molecular gas, h = 100 pc and the extincting 
layer is horizontally isotropic. The true extinction toward the inner Galaxy 
is most likely dominated by the molecular ring and nuclear region given the 
molecular gas distribution. However, precise estimate of the true distribution is 
not available and an horizontally isotropic model will adequately represent its 
systematic effect on the photometric distances. 

Of the more than 158,000 good flux-quality sources listed in IRAS PSG, 
5,736 satisfy both flux limit and variability criteria. Their spatial distribution is 
shown in Figure |l|. To obtain variability data, at least two epochs are needed. 
Unfortunately, IRAS’ multiple epochs did not have complete sky coverage. Most 
of the coverage (77% in the galactic plane) was achieved in HGON 2 and HGON 
3 separated by roughly 7.5 months on average. The rest of the galactic plane is 
poorly sampled (shaded regions in Figure |^). For this analysis, all the data in 
the poorly sampled sectors have been excised, reducing the size of the sample to 
5,500 stars. 


3. Method overview 

All of the selection effects but especially data incompleteness greatly 
complicate the analysis. Bayesian techniques are ideally suited to parameter 
estimation over data with general but well-dehned selection criteria and underlies 
both the maximum entropy and maximum likelihood procedures. Below, we 
will parameterize the source density by an exponentiated orthogonal series with 
unknown coefficients Aij and Bij (cf. eq. 0). In this context, the basic theorem 
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of the theory reads: 




P {{A,}, {B,,} \I)-P{D I {A,}, {p.,}, I) 
P{D\I) 


( 2 ) 


The probability P{{Aij}, {Bij} \D, I) is the conditional (or posterior) 
probability of the coefficients of the source density provided the data 
(D) and information (J) describing its incompleteness. The probability 
P ({Ajj}, {Bij} I /) is the prior probability (or simply, prior) of the coefficients 
provided only the information. Following Bretthorst (1990), we assign the prior 
using the maximum entropy principle. In our case it is constant implying that all 
coefficient values are equally likely initially. The function P {D \ {Aij}, {Pp}, I) 
is the direct probability which describes the likelihood of data given the 
coefficients. Finally, P {D\ I) is a normalization constant which may be omitted 
provided that the posterior probability is normalized. 

With these dehnitions, it follows that 


P {P,,} I P, /) = Const ■ P (P I {A,^}, {P,,}, I), (3) 

or in words, the posterior probability is proportional to the likelihood function. 
Therefore, the best estimate of posterior probability is obtained for the set 
coefficients which maximize the likelihood function. 

4. Likelihood function 

The likelihood is the joint probability of the observed stars given a source 
density. We may then consider the probability of observing a star with intrinsic 
luminosity in the range (L, L + dL) to be detected in the distance interval 
(s, s + ds), in the azimuth interval {1,1 + dl), in the galactic latitude interval 
{b, b + db) and with magnitude in the range (m, m + dm). Assuming a normal 



distribution of intrinsic luminosities L and a normal error distribution for the 
apparent magnitudes m this becomes: 

Pn (s, I, b, m, L\ am, ctl, K 12 , h, Rq) ds cos b db dl dL dm = 

0, ds cosbdbdl dL dm. (4) 

Here s, I, b are coordinates about the observer’s position, r, 0, z are coordinates 
about the center of the Galaxy, C is the normalization constant, S(r, 0, z) is 
the source density at galactocentric radius Rq, L and aL are the mean intrinsic 
luminosity and the dispersion of the sample, am is the measurement error in 
magnitudes and m = m{s,b) is given by equation (^). Alternatively, we may 
replace luminosity by absolute magnitude: 

Pn (s, I, b, m, M \ am, ctm, K 12 , h, Rq) ds cosbdbdl dMdm = 

G ■ S(r, 0 , z) ^-im-rnf/2al ^2 ^ ^ 5 ^ 

where M and ctm correspond to L and a^. The Gaussian distributions in L or 
M in the above two equations can be generalized to an arbitrary luminosity 
function for traditional star count applications. Although we will not give the 
general expressions below, the development is parallel. 

Since the convolution of two Gaussians is a new Gaussian whose variance is 
the sum of the two individual variances 

^m,eff A ; (®) 

equation (|) can be rewritten as 

Pn (s, I, b, m\ am,eff, k, H, Ro) ds cos b db dl dm = 

G ■ S(r, 0 , z) ^ 
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after integrating over the unmeasured absolute magnitude M. For notational 
clarity, we will omit the subscript “eff” and write simply am- The constant C is 
determined from the normalization condition: 

/ + 00 , _ 9 r rSmax(b) „ 

g-(m-m) / J dl J ds J S(r, 0, z) cosbdb = 1. (8) 

The integration over / runs over entire circle except missing azimuthal sectors, 
explicitly accounting for missing data at particular ranges in azimuth. The 
limiting distance Smax in the /, b direction incorporates the 2 Jy flux limit. 

In a standard star count analysis no explicit distance information is provided 
and s is eliminated from analysis by integration, yielding 

Pn {I, b, m\ ...) cos b db dl dm = 

fSmax{b) _ 2 /^^ 9 

C S(r, 0, z) ds cos b db dl dm. (9) 

Jo 

For our relatively small sample of IRAS stars, sensitivity to vertical structure 
will be poor. This motivates replacing the general unknown three-dimensional 
disk density with a density which depends on radial position and azimuth alone: 
S(r, 0, z) = E{r, 0). 

Finally, the joint probability of observing N stars selected from the IRAS 
PSC is 

N 

L = Ptotal = n ^riil, &, m| . . .). (10) 

n=l 

Expressing the likelihood function in logarithmic form, our desired solution is 
the set of parameters which maximize 

N 

logT = ^logPn(/,&,m| ...). (11) 

n=l 

This and nearly all star count analyses reduce to standard problem of 
density estimation: End the density function /(x), which satisfies non-negativity 
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constraint 

fix) > 0 (12) 


and integral constraint 

J f{x)dx = 1 (13) 

which best describes the observed data distribntion. Both parametric and 
non-parametric estimation techniqnes have been nsed to solve this problem (e.g. 
Silverman 1986; Izenman 1991). For inhomogeneons mnltidimensional data, 
the positivity constraint is cnmbersome. However, searching for the nnknown 
fnnction f{x) in the form of an exponentiated orthogonal series (Clntton-Brock 
1990), gnarantees positivity. A candidate stellar snrface density is: 

f '^max Jmax >. 

S(r, (j)) = exp EE [Aij cosjcj) + Bij sinjcj)] Jjikjr) L (14) 

i=l j=0 ^ 

where Jj(x) is Bessel fnnction of order and kj is root of Bessel fnnction 
til 

of order and are chosen to prodnce a complete orthogonal set over the disk 
of radins Rmax- The coefficients Aij,Bij are the parameters to be determined. 
There is no loss of generality in taking the Fonrier-Bessel series althongh the 
choice is arbitrary. 


5. Results 

5.1. Sensitivity to incompleteness 

A major advantage of the approach presented here over that in Paper I is 
that the significance of inferred strnctnre is robnstly qnantified. In particnlar, 
we can test the sensitivity of selection effects to the detection of a bar. To test 
the presence of the coverage gaps, we generated fonr sample disks of 1,000 stars 
each nsing the sonrce density (11^) with ^A‘fj + Bfj = 1 for j = 0, 2 and zero 
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otherwise and the following bar position angles: 0°, ±45°, and 90°. The root 
sum square of the coefficients Aij and Bij represents the strength of radial 
component for the polar harmonic. Figure shows the restored strength of a 
harmonic Afj ± Bfj as a function of the position angle of the bar. Insensitivity 
of these strengths to bar position angle suggests that missing azimuths will not 
obscure the inference of true bar. The computed values are consistent with the 
expected value of unity. 

Conversely, regions of missing data can produce non-axisymmetric 
distortions, and in principle, suggest the existence of a bar in initially 
axisymmetric sample. However, analysis of a simulated axisymmetric disk 
(Hio = A 20 = 1; all others = 0) and the same azimuthal incompleteness as in 
the real sample shows that the power in the non-axisymmetric harmonics is 
about 3% of the axisymmetric contribution. Together these tests suggest that 
the misidentihcation of a bar relative due to missing azimuthal sectors alone is 
unlikely. 


5.2. Application to IRAS data 

The formalism developed in §|| requires the distance to galactic center Rq, 
extinction in the plane K 12 and average luminosity of the AGB stars L. We 
adopted Rq = 8.0 kpc, Ku = 0.18 mag/kpc and L = ShOOL©. The method can 
be straightforwardly modihed for complex models (e.g. patchy or non-uniform 
extinction), the only limitation here is the CPU available and sufficient data to 
attain a satisfactory measure of conhdence. 


Choosing the truncation of the series in equation (|T^ poses a problem 


common to many non-parametric density estimations: because too few terms 
result in large bias and too many terms increase variance, imax, jmax would 
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be best determined by jointly minimizing the bias and the variance. However, 
this approach is computationally prohibitive due to the integral in equation 
@ and the normalization (H). Therefore, a heuristic approach was adapted 
in selecting imaxi jmax based on the increase in the likelihood function when 
a particular term or set of terms is added. Signihcance could be quantihed 
in terms of the likelihood ratio (Wilks 1962) but we have not done this here. 

In addition, the hardware available to us makes it impossible to sample the 
parameter space beyond imax = 4, jmax = 4. Nevertheless, up to that limit, the 
space was sampled thoroughly, with some of the solutions shown in Figure ^ 
along with the corresponding offsets of the likelihood function (the lowest value 
of likelihood is set to 0 for ease in comparison). Some of the hgures feature the 
ghost peaks due to the absence of data beyond the galactic center or in missing 
azimuthal sectors (see Figs. |l| and H). The likelihood analysis may attempt to 
place a non-existing source density peak in that region, provided it will increase 
the overall score. We will pursue penalizing the likelihood function and other 
procedures for choosing an alternative prior (dropping the assumption that all 
coefficients in (^) are equally likely initially) in future work. 

More importantly, all reconstructions in Figure ^ imply a jet-like feature 
in the hrst quadrant. As in Paper I, the depth of our sample (estimated to 
correspond to a mean distance of 7 kpc in the plane) prevents ascertaining 
whether this feature corresponds to a bisymmetric bar or is a lopsided distortion. 
However, decreasing the flux limit to 1 Jy leads to detection of similar feature on 
the far side of the Galaxy, suggesting a real bar. This motivates a reconstruction 
with enforced bisymmetry, shown in Figure Here the corresponding prior 
assigns zero values to coefficients of odd azimuthal order. The likelihood value 
(the origin is the same as in Figure ^) has dropped substantially, because the 
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resulting density lacks data support beyond the Galactic center. In both hgures, 
the bar is well dehned and has a similar length and position angle. 

To quantify the strength and position angle of the bar, we htted the 
isodensity contours {imax = jmax = 4) by ellipses. The logarithm of a suitable 
likelihood function for estimating the semi-major axes, eccentricity and position 
angle is 


logT = ^ 


Mr -|2 


(15) 


2 = 1 ^ 

where 0) is the reconstructed density function and C is isodensity level. 

The summation runs over equally spaced points on ellipse. For a given ellipse, a 
grid of semimajor axis values are specihed and the surface density C, position 
angle 0o and eccentricity e which maximizes logL are found. The results are 
presented in Figures ^ and |^. 

Figure || indicates that the density prohle drops to half of its central value 
at about 4 kpc. The half-length would then be about 4 kpc, in good agreement 
with the value obtained in Paper I. If we take this value as the size of the 
major axis of the bar, then the axis ratio varies from 2.2 in the central regions 
to 2.7 in the outer regions of the bar. The value of the position angle for 
the entire extent of the bar (out to 4 kpc) is ~ 19°. The accuracy of the 
position angle determination can be quantihed in terms of conhdence interval, 
making use of the fact that in the limit of large number of sources N, the 
likelihood in n dimensions is distributed as x^/2 with n degrees of freedom (e.g. 
Lehmann 1959). We analyzed the likelihood as the function of a single variable 
~ orientation angle of the bar in the plane. The analysis gives the uncertainty of 
1° at 3cr level. 

Another way to determine the parameters of the bar is to look at the map of 
the ratio of non-axisymmetric to axisymmetric components of the density. The 
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ratio displays two peaks at 3.3 ±0.1 kpc located on the opposite sides from the 
center, the line connecting them has the position angle of ~ 24° ± 2°. The peak 
ratio, the relative strength of the bar, is 0.73. This implies the existence of a 
strong bar in the intermediate age population responsible for the AGB stars. 

5.3. Disk scale length 

Having calculated the source density, we are in a position to characterize the 
parent population of the IRAS variables. In Paper I, we assumed that these 
variables represented a disk population based on their flux distribution but 
several colleagues have suggested in discussion that the IRAS variables are more 
likely to be bulge stars. Here, we determine the scale length of the population 
in the Galactic plane. For comparison, we £t our reconstruction by an oblate 
spheroid model (GO bulge model from the DIRBE study by Dwek et ah 1995): 

SGo(x,|/) = Soe-°'5^', (16) 

with = (x^ ± ip‘)/rQ. The scale length tq is found by minimizing the following 
cost function while simultaneously satisfying the overall normalization constraint 
for Ego (eq. [T|): 

cost = J (fr 

To estimate the value of rg, we used the covariance matrix from the likelihood 
analysis used to determine E^ec to make 5000 Monte Garlo realizations of the 
source density. The ensemble of realizations, then, have E^ec as their mean. 
For each realization, we found vq by minimizing the cost function (|I7D and 
the resulting distribution of scale lengths is shown in Figure Our result 
rg = 4.00 ± 0.55 kpc indicates that the IRAS variables have the scale length of 
the old disk population. This value is in good agreement with the scale length 4.5 
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kpc reported by Habing (1988), derived from analysis of a color-selected IRAS 
sample. Dwek’s value obtained by analyzing bulge emission was tq = 0.91 ± 0.01 
kpc. The factor of 4 difference between the scale lengths suggests that the IRAS 
bar and the bulge-bar belong to distinct populations. 


5.4. Optical depth due to microlensing 


Originally proposed as a test for dark matter in the Milky Way halo 
(Paczyhski 1986), gravitational microlensing was later shown (Griest et al. 1991; 
Paczyhski 1991) to be potentially useful for extracting information about the 
inner regions of our Galaxy. Three groups (OGLE, MAGHO and EROS) are 
monitoring stars in the Galactic bulge for gravitational microlensing and have 
found higher event rates than most theoretical estimates. Udalski et al. (1994) 
derived lensing optical depth r = (3.3 ± 1.2) x 10“® toward the Baade’s window 
(/ = l°,b = —3.9°) based on the analysis of the OGLE data, and MAGHO group 
reported r = 3.9lJ;2 x 10“® (Alcock et al. 1995a) estimated from the sample 
of clump giants, while theoretical estimates give optical depths in the range 
0.5 — 2.0 X 10“® (e.g. Alcock et al. 1995a; Evans 1994). Following Paczyhski’s 
et al. (1994) suggestion that a bar with a small inclination angle could enhance 
the optical depth, Zhao et al. (1995) have developed a detailed bar model and 
found r = (2.2 ± 0.5) x 10“®. Here, we estimate the optical depth using our 
density reconstruction, S^ec, assuming that our AGB sample represents the 
entire stellar disk. 

The lensing optical depth is dehned as the probability of any of the sources 
being lensed with magnification factor A > 1.34, with 


A 


u^ + 2 


r 



u^/u^ -|- 4’ 


(18) 
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(Refsdal 1964), where r is the distance between the projected position of the 
sonrce and the lensing mass, Re is the radius of Einstein ring. Kiraga & 
Paczyhski (1994) derived 


r = 


4iG r f/o'’- P doA p Dl*-‘P dD, 


(19) 


c2 

where Dg is the distance to the sources, is the distance to the deflectors 
and the free parameter (3 accounts for detectablity of sources in a flux-limited 
survey. The reasonable range is —3 < /3 < —1 and we take (3 = —1 following 
Evans (1994) and Kiraga & Paczyhski (1994). The density p = pbuige + Pdisk-, 
with Pbuige given by equation (1) of Kent (1992), and 




( 20 ) 


where S 44 is the surface density of our i = 4,j = 4 model (p!4D and h = 0.325 
kpc is the scale height. We explored two possible normalization prescriptions: 
(1) Assign a local column density of ~ 50 MqPc“^ (“canonical disk” following 
Kuijken & Gilmore 1989; Gould 1990). The mass of the disk in this case is 
Mdisk = 1-95 X IO^^Mq. (2) Assign the total disk mass of M = 6 x 10 ^°Mq 
(B ahcall & Soneira 1980). The second normalization gives local column density 
of approximately 100 MqPc“^ (“maximal disk” of Alcock et al. 1995b). We 
prefer the latter here because the optical depth estimate depends on the global 
mass distribution rather than the local density. In addition, there are some 
indications that the variation of the column density with galactic longitude may 
be quite significant - a factor of 2 — 3 (Rix & Zaritsky 1995; Gnedin, Goodman 
& Frei 1995). The mass of the bulge is Mbuige = 1-65 x IO^^Mq. 

For the canonical disk case, the total lensing optical depth at Baade’s 
window is 1.1 x 10“®, and both bulge and disk lenses contribute 50% to that 
number. Most of the optical depth (76%) is due to lensing of bulge sources. If 
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the disk is maximal, optical depth is 1.6 x 10“®. Disk lenses now account for 
1.1 X 10“® (68% of the total optical depth) and the contribution by bulge sources 
still dominates (59%). For both scenarios, optical depth is a function of the 
orientation of the bar. We investigate the enhancement produced by the bar over 
axisymmetric models of the disk p cx where R = 3.5 kpc for hxed 

disk mass. Figure displays the ratio of optical depths of non-axisymmetric 
to axisymmetric disk models as a function of the position angle of the bar for 
both normalization scenarios. The difference between the two curves illustrates 
the role of the disk in lensing. The largest enhancement of approximately 30% 
obtains when the bar is aligned along the line of sight as expected. The ratio of 
optical depths decreases gradually when the bar is in the hrst Galactic quadrant, 
with > 20% enhancement out to 0o = 50°. 

Current generation optical-band lensing surveys have concentrated on 
low-extinction bulge-centered windows to maximize the lensing event rate. 

An infrared-band lensing microlensing survey would be less constrained by 
extinction and therefore more efficient probe of the overall structure of the 
Galaxy. In particular, any bar which is not perfectly aligned along the 
Sun-Galactic Center axis will produce an asymmetry in the optical depth. 

We describe this asymmetry by the ratio of the difference in optical depths at 
positive and negative longitude to their arithmetic mean. This ratio is shown in 
Figure 0 for our model (cf. eqns. ^ and ^). Comparison with the Bahcall 
& Soneira model (1980) suggests that /d ~ —1 is a fair approximation of the 


high-luminosity end of the disk luminosity function. Therefore, equation (|l^) 


also applies at large |/| where both lenses and sources are disk members. The 
large 40% asymmetry about |/| 30° is due to a local increase in the surface 

density at negative longitudes close to the observer (Figure More important 
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than the details of asymmetry is the suggestion that a pencil-beam microlensing 
survey in the infrared would be sensitive to global asymmetries in the stellar 
disk component. Confusion is not a limitation at 6 = 0° for larger values of |/| 
and the optical depth has a magnitude similar to Baade’s window. 

6. Summary and discussion 

This paper explores a model-independent Bayesian estimation of the stellar 
density from star counts, rigorously accounting for incomplete data. The 
general approach can incorporate multiple colors and even different databases. 
The usual high dimensionality and topological complexity of the posterior 
distribution, however, complicates both optimization algorithms and subsequent 
moment analyses. We propose here a hybrid downhill plus directed-search 
Monte Carlo algorithm; the former speeds convergence and the latter facilitates 
the location of the global extremum. Other similar and potentially more efficient 
techniques which can bypass the extremization step altogether (such as general 
Markov Chain Monte Carlo) are worth careful consideration. 

Application of the technique to the variability-selected sample described 
in Weinberg (1992), assumed to be AGB stars, conhrms the presence of a 
strong non-axisymmetric feature in the hrst Galactic quadrant. By imposing 
bisymmetry on the source density, clear signature of a bar is obtained. The size 
and shape of density isophotes suggests a bar semi-major axis of approximately 4 
kpc and position angle of 0o = 18° ±2° at the outer edge of the bar. The analysis 
of the scale length for the AGB candidate distribution gives ro = 4.00 ± 0.55 
kpc, indicating that these objects are part of the old disk population. 

Finally, we use our estimate for non-axisymmetric Galactic disk to explore 
the dependence of optical depth to gravitational microlensing by bulge and disk 
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stars. The disk bar does enhance the optical depth r towards Baade’s window 
by roughly 30% but the overall value is still roughly a factor of two below the 
MACHO result r = x 10“®. Of interest for future microlensing surveys is 

the hnding that our inferred large-scale bar will produce a signihcant asymmetry 
in r at positive and negative longitudes beyond the bulge. The peak asymmetry 
for our model occurs at |/| = 30° and at 6 = 0 we predict similar values of r 
to the Baade’s window held. Such a survey might best be carried out in the 
infrared to take advantage of the low interstellar extinction and colors of the 
late-type giants. At |/| 30°, confusion should not be a limitation at 6 = 0°. 

We thank Steve Price and Mike Skrutskie for comments. This work was 
supported in part by NASA grant NAG 5-1999 and the Alfred P. Sloan 
Foundation. 


A. Luminosities of AGB stars 

The luminosities of AGB variables and the inference of their progenitor 
masses plays a role in constraining the stellar evolution history of the Galaxy 
and has received some attention. Investigations based on theoretical approach 
(Iben & Renzini 1983) and observations of sources close to the Galactic center 
(Jones & Hyland 1986) placed the luminosities somewhere between a few 
xlO^L© and 6 x lO'^L©. Van der Veen & Habing (1990) revised the results 
of Jones & Hyland based on the analysis of a larger sample of OH/IR stars 
and found that the luminosities are in the range 10^ - with the peak of 

the distribution about 5, 000 — 5, SOOT©. They suggested the variability of the 
sources (Am < 2*”) and possible selection effects as main reasons for higher 
limits of Jones & Hyland. They also noted that as many as 20% of the stars 
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may be in the low-luminosity tail of the distribution but only 2% or fewer can 
exceed the upper limit. Kastner et ah (1993) obtained kinematic luminosities 
based on the radial velocities of circumstellar envelopes with respect to the 
LSR and distances derived from the Galactic rotation curve. They found the 
range of 1.3 x 10^ - 2 x IO^Lq with average uncertainty of factor of 2. The 
theoretical estimate was recently revised by Groenewegen et al. (1995) who 
obtained luminosity functions for carbon and oxygen-rich stars based on the 
synthetic evolution. They found a mean luminosity for Galactic carbon and 
oxygen-rich AGB stars to be TOSOL© and 3450Lq, respectively. They stated 
“the luminosity of a typical Galactic AGB star is in any case less than the IO^Lq 
often assumed”. Habing (1988) reported the average luminosity of 4000Lq for 
a color selected sample from IRAS PSG catalog. Finally, analysis of a sample 
of oxygen Miras using P-L relation established on the observations of LMG 
Miras (Feast et al. 1989), places their average luminosity at L = 3900 ± 450Lq. 
Unfortunately, we can not use the P-L relation, since IRAS had insufficient 
temporal coverage to reliably constrain periods. Rather, we approximate the 
source density by an axisymmetric distribution at i?o = 8 kpc and choose the 
average luminosity which maximizes the likelihood function. The results for 
different number of radial terms are shown in Figure |T^. For ten terms, the 
maximum likelihood of this axisymmetric density is achieved when L k. OOOOLq. 
We adopt L = OSOOLq which is the low end of published results and interpret 
our statistical analysis as a consistency check. 

B. Computational Notes 

Likelihood maximization is the rate limiting step in inferring the surface 
density from a source catalog. The cost of computing the likelihood is 


proportional to the sample size so analyses of very large data sets will be 
technically challenging. Our “workhorse” algorithm for locating the maximum 
of the likelihood function is the conjugate gradient method which is thoroughly 
discussed in the literature (e.g. Press et ah 1988). We have adopted an 
implementation by Shanno & Phua (1976, CONMIN). The algorithm has good 
convergence properties, but requires a good initial approximation. Near the 
expected quadratic maximum the convergence should be extremely rapid. 

However, the likelihood function may have a large number of extrema, 
limiting the use of the standard downhill technique. In such cases, the Simulated 
Annealing (SA) algorithm (Metropolis et ah 1953; Otten & van Ginneken 1989) 
has the advantage. It places no restrictions on continuity and easily incorporates 
arbitrary boundary conditions and constraints. Adaptive Simulated Annealing 
(ASA, Ingber 1989)—a faster version of the SA algorithm—proved to be effective 
in narrowing the domain of the search to the comparatively small region in 
parameter space. However, in the vicinity of the extremum it converges slowly. 

The complementary features of the two techniques, suggest the following 
two-step hybrid scheme: 

1. Use a directed search algorithm (ASA) to isolate the global maximum. 
Although SA class of algorithms converge slowly, there is a probabilistic 
guarantee of convergence: the probability of hnding the maximum is 
inversely proportional to the total number of iterations to some power (e.g. 
Shu & Hartley 1987; Ingber 1993). 

2. After either a limiting number of steps or a signihcant drop off in 
convergence, use the current ASA solution as input to conjugate gradient 
scheme. This is motivated by our expectation that the true maximum of 
the likelihood function will be a quadratic form in the unknown variables. 
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This sequence can be repeated again, in case if Step 2 fails to find a well-defined 
maximum. The scheme is difficult to analyze but appears to work well in 
practice and is potentially useful for large parameter space and complex 
geometry (boundary conditions, irregular likelihood function) cases. 

The entire computation time scales as the number of coefficients M (total 
number of Aij and Bij in the sum in eq. |^) and the sample size N: N{2M + 1). 
Computation of the Hessian matrix requires CPU time proportional to M^N. 
For a large M, this is the bottleneck. However, the algorithm is straightforwardly 
parallelized by partitioning the data. 
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Fig. 1.— The sample of 5,500 IRAS PSC variables (dots). The Sun is located 
at X = —8, Y = 0. The data from the shaded sectors are eliminated from 
the analysis. The circle shows the distance in the plane where an AGB star 
(L = 35 OOL 0 ) can be detected. 

Fig. 2.— The same sample projected on the X-Z plane. All the data are inside 
the region bounded by two solid lines which are solutions of the equation (^j. 

Fig. 3.— The amplitude of harmonic coefficients as functions of the position 
angle of the bar. Open triangles: i = 1, j = 0; open squares: i = 1, j = 2] hlled 
triangles: i = 2,j = 0; hlled squares: i = 2,j = 2. The symbols are slightly 
offset along the x-axis for clarity. 

Fig. 4.— The reconstructed density prohles. Ten equally spaced contours 
between 10% and 100% of peak value are shown in each panel. 

Fig. 5.— The reconstructed density prohle obtained with assumption of 
bisymmetric source density. There are 10 contours between 10% and 100% of 
peak value in each panel. 

Fig. 6.— Isophotal hts to the reconstructed source density: surface density C 
normalized to its central value (left scale, solid line) and axis ratio a : b (right 
scale, dashed line) versus semimajor axis. 

Fig. 7.— The position angle cf) in degrees (left scale, solid line) and eccentricity 
of ellipses (right scale, dashed line) versus semimajor axis. 

Fig. 8.— The distribution of the scale lengths ro in 5000 realizations of the source 
density (histogram). The best ht normal distribution is shown (solid curve) with 
mean and rms value as labeled. 
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Fig. 9.— The ratio of optical depths toward Baade’s window obtained with non- 
axisymmetric (bar) and axisymmetric disk models as the function of the position 
angle of the bar, 0o- Solid line - “maximal disk”, dashed line - “canonical disk” 
(see text). 

Fig. 10.— Asymmetry in the microlensing optical depth. The disk is “maximal”. 
Solid line, dashed line and dotted line represent cuts with b = 0°, 2° and 4°, 
corresp ondingly. 

Fig. 11.— Average optical depth as the function of the galactic longitude. The 
lines represent the same latitudes as in Fig. 

Fig. 12.— The values of the likelihood function with varying luminosity of the 
sources. The center of the Galaxy is hxed at 8 kpc. Solid line — imax = 2, dotted 
line — imax = 6, dashed line — imax = 10. 




































































